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Abstract. We consider stochastic rules of mass transport which lead to steady states 
that factorize over the links of a one-dimensional ring. Based on the knowledge of the 
steady states, we derive the onset of a phase transition from a liquid to a condensed 
phase that is characterized by the existence of a condensate. For various types of weight 
functions which enter the hopping rates, we determine the shape of the condensate, its 
scaling with the system size, and the single-site mass distribution as characteristic static 
properties. As it turns out, the condensate's shape and its scaling are not universal, but 
depend on the competition between local and ultralocal interactions. So we can tune 
the shape from a delta-like envelope to a parabolic-like or a rectangular one. While 
we treat the liquid phase in the grand-canonical formalism, we develop a different 
analytical approach for the condensed phase. Its predictions are well confirmed by 
numerical simulations. Possible extensions to higher dimensions are indicated. 
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1. Introduction 

Stochastic mass transport describes the movement of a generic mass from one space point 
to another as a classical stochastic process. "Mass" stands for a generic quantity that 
is conserved throughout the process. It can be as different as a vehicle in macroscopic 
traffic or a macromolecule moving along the cytoskeleton pQ. Above a certain density, 
mass condensation may be observed, so that even in the limit of infinitely many masses 
a finite fraction of all masses piles up at a certain site, although the dynamic rules are 
fully symmetric. This is an expression of spontaneous symmetry breaking. In these 
systems it can happen even in one dimension, as they are in general out-of-equilibrium. 
When the mass transport describes vehicles, mass condensation corresponds to traffic 
jams. Other applications of condensation are found in granular flow j2], clustering [3], 
they are manifest as gelation in networks where a single node takes a finite fraction of 
all links [4], as Bose-Einstein condensation, or in phase transitions of quantum gravity 

Many of these systems can be modeled as a set of particles occupying discrete levels, 
or "boxes" on a one-dimensional grid. The balls- in-boxes model (B-in-B) [6], its non- 
equilibrium version, the zero-range process (ZRP) [7], or a more general model [8] with 
continuous masses are well-known examples. Although these models are more abstract 
than realistic ones, they serve as paradigms since the stationary state can be derived 
analytically as fully factorizing over the sites of the grid, or, more generally, over the 
nodes of an arbitrary graph [9]. The factorization is due to ultra-local ("zero-range") 
rules assumed to govern the dynamics of particles. This makes the phase structure 
of the systems accessible and enables comparison with experiments on condensation 
phenomena, at least on a qualitative level. 

As a natural consequence of the lack of interactions between particles at neighboring 
sites, the condensates in the B-in-B or ZRP models — when they occur — they always 
occupy a single site. The question then arises as to how the shape of the condensate 
changes in the presence of interactions between sites that tend to flatten out the 
condensate's profile but still preserve translational symmetry and conserve the current. 
A first answer has been given in [10] for a model that is related to a solid-on-solid (SOS) 
model [TT1 [T2] , supplied with dynamical rules that drive the system out-of-equilibrium. 
There, the steady state factorizes over pairs of sites of a one-dimensional ring topology 
which allows for nearest-neighbour interactions while making the system analytically 
solvable. The condensate is extended over a range that scales as y/~N if N is the system 
size. The hopping rates were determined by weights g(m,n) that factorize into an 
ultralocal part, depending on the occupation numbers of a single site, and a local part, 
depending on the occupation number of neighboring sites. For this model, the scaling 
of the condensate's extension W with yN can be derived as the first-return time of a 
random walker after W time steps (see [10] ). This is possible when the very formation 
of the condensate is described in terms of a random walk, for which the walker chooses 
his step sizes according to a distribution that is given by the weights g(m, n). One may 
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naively expect that this analogy carries over also to other choices of weight functions, 
so that an extension with \/N would be a universal feature of condensates with PFSS. 

In this paper we shall show that this is not the case, i.e., the scaling is non-universal. 
Instead we can actually tune the shape of the condensate and its extension via different 
appropriately chosen weight functions. In terms of a random walk, the walker would no 
longer use the weights g(m, n) as distribution that determines his choice of step sizes, so 
that we have to look for another derivation. As it turns out, the relevant characteristics 
for both local and ultralocal weight factors is the range of interactions in occupation- 
number space: we call short-range interactions local weight factors K(x) = K(\m — n\) 
and ultralocal factors p(m) which decay faster than any power in their argument, 
and long-range interactions those which decay as a power of their argument (x or 
m, respectively). For local interactions with exponential decay (short-range) and an 
ultralocal part that approaches a constant above some occupation number m max , we 
are able to analytically derive the static characteristics of the condensate: the critical 
density for condensation, the extension of the condensate for a given size, the shape of 
the condensate, that is its average occupation number as a function of spatial position, 
the fluctuations of the condensate and therefore also the mass distribution, that is the 
probability for finding m unit masses at a given site of the ring. The derivation is not 
restricted to the fluid phase for which the description in terms of a grand-canonical 
partition function holds up to the critical mass density. In our approach the condensed 
phase becomes analytically accessible due to the assumption that the partition function 
factorizes into a part Z c that describes the contribution of the condensate, and a critical 
background Z b , for which the uniform particle distribution has just the critical density, 
so that all excess mass is absorbed by the condensate. The factor Z c can be explicitly 
calculated when it is independent of the ultralocal part of the weight, i.e., it is of 
p(m). This happens in particular for the choice of weight functions as it was used 
in [10]. Otherwise we use what we call the fixed-envelope approximation. Here we 
approximate Z c via the probability P(W) of having a condensate extended over W 
sites with unknown but fixed envelope h(t) for which we use a form that is inspired by 
numerical simulations of the condensation. The condensate's extension follows in these 
cases from the maximum of the probability P(W). 

As it turns out, within the same fixed-envelope approximation scheme, when local 
and ultralocal weights are "short-ranged", we can tune the extension to scale with the 
system size as N a with < a < 1/2. When both parts of the weights are long-ranged, 
the condensate gets localized to a single site as in the ZRP. For long-range local part and 
short-range ultralocal part, the condensate takes a rectangular form, its height scales 
proportionally to the system size, while its extension remains constant: features that 
remind to finite-size scaling of a first-order phase transition. Therefore hopping rates, 
leading to PFSS, do not necessarily lead to extended condensates, and if they do so, the 
shape of the condensate is non-universal. 

The paper is organized as follows: In the next section we introduce the model 
and derive its stationary state from the corresponding master equation. In section 



Mass condensation in one dimension with pair-factorized steady states 



4 



|3J we discuss short-range interactions with an ultralocal part that is constant above 
some threshold m max . Under mild assumptions on the factorization of the partition 
function we derive the critical density for condensation, the scaling of the extension of the 
condensate, the condensate's shape and its fluctuations, and the single-site distribution 
of particles, in particular also in the condensed phase. In sections H] and [5] we use 
a more heuristic argumentation to predict the shape of the condensate for all other 
choices of weight functions, in particular for power-like decay of local weights. Section 
[6] contains the summary and conclusions. Some supplementary material is provided in 
the Appendices. 



2. The model and its stationary states 

We assume that we have M particles placed initially at random on N sites forming a 
closed chain. Each site % carries a certain number rrii of particles, which may range 
from zero to M. The periodicity implies that m^ + i = m\. With rate tt(mj|mj_i, 
the particle jumps out of site i it occupies, and moves to one of the two neighboring 
sites. The hopping may be asymmetric: with probability r the particle departures 
to the right, and with probability 1 — r to the left. In a computer simulation, this 
corresponds to choosing at each time step the departure site at random and, if it is not 
empty, moving a single particle with probability proportional to w(mj|mj_i, m i+1 ) (after 
a proper normalization) to one of the neighboring sites. Furthermore, we assume that 

u{mi\mi-x, m i+1 ) = /(m, i; m i _ 1 )f(m i) m i+1 ), (1) 

with some non-negative function f(m,n). The hopping rate is thus a product of two 
factors, one for each nearest neighbour. As we shall see, this is the crucial assumption 
that leads to the pair-factorized steady state. The model is very similar to the one from 
[TO] , the only difference is that there the hopping rate was a more general product of 
two factors /i(mj,mj_i) and f2(n r H, ^i+i)- Some aspects of the model were analyzed 
in [15] . In what follows we shall discuss this choice more thoroughly and study some 
special cases. 

At each time t, the state of the system is specified by the set of occupation numbers 
rh = {mi,...,mjv}. Let us denote by P(rh,t) the probability of having the system in 
a particular state fh at time t. The evolution of P(fh,t) is governed by the master 
equation: 

dP (™^ = [W(rh' -> m)P(m', t) - W(m -> fh')P(m, t)] , (2) 

(J/L _ ( 
m! 

where W(m! — > fh) is the transition rate from state fh! to state fh. In this paper we are 
interested in the steady-state solution of that is in time- independent probabilities 
P(fh) = P(rh,t — ► oo). As follows from (T5]), 

^P(m)W(m -> fh') = ^P(m)W{m -> rh). (3) 
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If both sides were equal separately for each fh', this equation would simply correspond 
to the detailed balance condition and hence to the equilibrium case. We know, however, 
that for r ^ 1/2 the system cannot be at equilibrium, because there is a net current 
of particles to the right (r > 1/2) or to the left (r < 1/2). It turns out, however, that 
the solution of ([3]) does not depend on the hopping asymmetry r, and takes the same 
pair-factorized form either out of or in equilibrium: 

N 

P(fh) = Y\_g(rni,m i+1 )5 mi+ ... +mN ,M, (4) 

i=i 

provided that g(m,n) is some symmetric function g(m,n) = g(n,m) and 

f{m,n) = — . (5) 

g(m,n) 

In what follows, we shall skip the Kronecker delta function 8 mi+ ... +mNj M for brevity. 
To prove (jl]), let us observe that the transition probability W(fh — > fh!) is non 
zero and reads u{rrii\mi-\, mj+i) only if fh' = {. . . ,rrii-\ — 1, rrii + 1, . . .} or fh' = 
{. . . , mi + 1, rrii+x — 1, . . .}, i.e., there is one more particle in the configuration fh' on an 
arbitrary site i and one less on one of its neighbors. From ([3]) we obtain: 

} u u(rni\mi-i, m i+l )P(m) 

i 

= ^2 [™( m i + M m i+i ~ l,mi_i)P(. . .,m; + l,m i+1 - 1, . . .) 

i 

+ (1 - r)u{m i+l + l|m» - l,m f+2 )P(. . . ,m, - l,m i+ i + 1, . . .)] . (6) 
We proceed by inserting the guessed pair-factorized steady state (jl]), and 

/ i \ gijm - g(mi - i,m i+1 ) 

u(mi\mi- 1 ,m i+1 ) = r — , (7) 

g(mi,mi-i) g{mi,m i+1 ) 

as follows from ([T]) and (jSJ), into equation (JBJ). After canceling some terms we obtain 

^ g{m-i, m - l)g{rrii - 1, m i+1 )g(m i+1 , m i+2 )Ri (8) 

i 

= r^2 g( m i-u l mi)g(m i , m i+1 - l)g(m i+1 - 1, m i+2 )Ri (9) 

i 

+ (1 - r) 2j£c(mi_i,TOi - l)g{rrii - l,m i+1 )g(m i+1 ,m i+ 2)Ri, (10) 

i 

where we used the symmetry of g(m,n) = g(n,m) and Ri = Ylj^{i-i i i+i] g( m ji m j+i)- 
Both terms without the factor r cancel immediately. The remaining terms proportional 
to r cancel after shifting the index i — >• i + 1 in the last sum. We have therefore proved 
that when the hopping rate takes the particular form ([7]), the system possesses the 
steady state factorized over pairs of neighboring sites. The proof can be extended to an 
arbitrary graph [16J. 

From the above discussion one sees that the hopping rate may be defined by the 
two-point weight function g(m,n). It is this very weight function that we shall vary 
throughout the following sections. In general we shall factorize it into a local interaction 
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factor K(\m — n\), and an ultralocal factor p(m) (playing the role of an on-site potential) 
according to 

g(m,n) — K(\m — n\)\/p(m)p(n), (11) 

where both K(x) and p(m) are some positive functions of x and m, respectively. When 
K(x) = 1, g(m,n) factorizes and we recover the ZRP with the weight p(m) since every 
>Jp(m) appears twice in the product over sites in (pi}. As indicated in the introduction, 
we call it short-range interactions in occupation-number space when K(\m — n\) or p(m) 
decay exponentially, and long-range interactions when they decay with some power of 
the argument. 

The most intriguing phenomenon in the ZRP is certainly condensation when a finite 
fraction of all particles occupies a single site. Condensation is triggered by the on-site 
potential p(m), which acts only on particles being at the same site. When we consider 
in the following PFSS with weights K(x) depending on the difference of occupation 
numbers x = |mj — at neighboring sites i and i + 1, we can explore the effect 

of nearest-neighbour interactions on the very condensation. Therefore, in what follows 
we shall determine the critical density for condensation and the static properties of the 
condensate for various choices of K(x) and p(m). 



3. Short-range interactions 

In the first part of this section we will assume the following weights: 

K(x) = e~ Jx , p{m) = e u5m -°, (12) 
with parameters J and U, proposed in [10]. The steady state reads: 

P{m u . . .,m N ) = exp (-J ^ \rrii - m i+1 \ + U ^<5 mi , j . (13) 

The delta function, reflecting the conservation of particles, was dropped for brevity. This 
will be our flag example, although we discuss later a more general case for which K(x) 
is an arbitrary function that falls off faster than any power law, and p(m) is constant 
for m greater than some m max . From ( fi~3]) we see how to interpret J and U. Firstly, 
because \rrii — m i+ i\ measures the rate of change of with i, fast-changing profiles 
are suppressed the more, the larger J is. Thus, J gives a kind of surface stiffness, 
which tends to flatten out the profile. Secondly, the term £V 5 mij0 measures the total 
number of free sites. The larger U is, the more probable are configurations with many 
empty sites. An interesting observation follows from comparison to the solid-on-solid 
(SOS) model [TT] . In the SOS model without pinning potential, the probability of a 
microstate is essentially given by (TT3"]) with U = 0. Such a model does not have a phase 
transition in ID. We will see below that it is the insertion of the [/-term that leads to 
a liquid-condensed phase transition, even in the ID equilibrium system. 
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3.1. Critical density for condensation 

In order to study the properties of the steady state we can ignore the fact that the system 
is out-of-equilibrium, and treat it with conventional methods of statistical mechanics. 
Thinking of P(m) as a probability of a microstate, we define the canonical partition 
function: 

Z(N,M) = ^IJ^^i)^^ (14) 

{rrn} i 

as well as the grand-canonical one 

Z N (z) = $>(iV, M ) zM = E ^ Jlgfatm*!), (15) 

M {rrn} i 

which is just the discrete analog of a Laplace transform of Z(N, M), and the fugacity z 
is determined from 

p=N\2: m *) = N^z—- (16) 

The partition function Zn(z) grows monotonously with z, and so does its derivative. 
Thus the left-hand side of ( fl6l) grows also with z. If the radius of convergence of 
Zn(z) is infinite, then for any finite p there exists some z > which obeys (1161) . Both 
ensembles, the canonical and the grand-canonical one, are then in our case equivalent in 
the thermodynamic limit. We will use this fact to calculate the distribution of particles 
7r(m): 

that is the probability of having m particles at a randomly chosen site. Taking the 
average in the grand-canonical ensemble we obtain 



7i(m) . . / , T mm2 T m2rn3 ■ ■ ■ Tm N mi (18) 

x * m2,...,mjv 

Z.N (■2') ^ r ^m\rri'2~^mimz ' ' ' T-vn,j^m\ IV 'l'\z) , (19) 



mi,...,mjv 



where 



T mn = z^ +n ^ 2 g(m,n). (20) 
If we define now m to be a normalized eigenvector of T mn to the largest eigenvalue 

Amax) 

^ Tmn&n A max 0m,; (21) 

we obtain for large iV that Z N (z) = A^ ax and ir(m) = 4> m . The eigenvector <p m has to 
decay with m. Otherwise, p calculated from ( fl6l) would be infinite in the thermodynamic 
limit. By analogy to the ZRP we can thus say that the system is in the liquid state - 
there is no condensation. 
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On the other hand, if Z^(z) has some finite radius of convergence z c , the derivative 
in ffl6l) can either grow to infinity for z — > z c , or tend to some constant. In the first 
case we have again no condensation, because for any p there is some real z < z c which 
obeys (|T6|) . In the second case, there exists a critical density 



with <f) m being now the eigenvector for z = z c , above which the grand-canonical ensemble 
does not exist. This in turn indicates a phase transition from the liquid to the condensed 
state. 

In order to calculate the critical density, it is convenient to assume that z c — 1. If 
z c 7^ 1, we can always shift it to one by rescaling the weight g(m, n). If <j) m and A max are 
now the eigenvector and the maximal eigenvalue of T mn (z = 1) = g(m,n), respectively, 
the recipe for obtaining the critical density is i) to find the eigenvector of g(m, n) to the 
largest eigenvalue, ii) to square its elements, iii) to find its mean value treating 0^ as 
probabilities. This can be done by numerical diagonalization, truncating the matrix to 
a finite size and keeping track of the size dependence. 

In some cases we can find the eigenvector analytically or, at least, we can decide 
whether p c is finite or not. First, let us state that because g(m,n) > for all m,n, 
from the Frobenius-Perron theorem all entries <p m must be non-negative. Moreover, 
the eigenvector to the largest eigenvalue is not degenerated. Since g(m,n) = g(n,m), 
all eigenvalues are real, and eigenvectors to different eigenvalues are orthogonal. This 
means that all other eigenvectors must have some entries negative, thus if we find a 
semi-positive vector that obeys the equation: 



with some constant A max , it has to be the one to the largest eigenvalue. 

If g(m, n) = g(M—m, M—n), it follows that m = 4>M-m and hence p c = M/2 oc iV 
tends to infinity for iV — > oc. This is the case of the SOS model, where g(m,n) = 
e -j\m-n\.^ there is no condensation, regardless of how big the density p is. We see the 
importance of the [/-term in (Tl3|) . Its role is to break the symmetry in g(m,n) under 
the change (m, n) — > (M — m, M — n). When the symmetry is explicitly broken in the 
occupation number space, from ( 1221) it follows that m must decay faster than ~ mT x 
if the critical density has to be finite. 

With the above remarks in mind, let us sketch the derivation of p c for the weight 
f fT2|) . First, z c = 1, so we do not need to shift the weights. As suggested by the functional 
form of g(m,n), we assume that the eigenvector m oc e M m,o+ Bm w ith some constants 
A,B. Indeed, by inserting it into ( 1231) with g(m,n) = e - J \^-n\+u(& m ,o+& n ,o)/2^ one g n( ^ s 
that the eigenvalue equation is fulfilled as long as 




(22) 



rn 




(23) 



n 



A = U/2, 



(24) 
(25) 
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One sees that when J < J , where Jq = U — ln(e u — 1), the eigenvector grows with m, 
therefore the critical density is infinite. This means that condensation is possible only 
for J > Jq. The unnormalized eigenvector reads <p m oc e u& ™.fi/i-{J-Ja)m , and hence the 
critical density becomes 

PC YZ=A (e^-e-^o))( e ^o)_i)- ^ > 

When J — > J , the above formula reduces to p c ~ l/(e 2 ( J_J °) — 1), the result obtained 
in [TU] where the asymptotic form <ft m ~ e Bm was assumed. As a byproduct we obtain 
the largest eigenvalue: 

A max = e u + e2J(1 _ e _ [7) _ 1 - (27) 
3. ,2. Extension of the condensate 

When the density p exceeds p c , the condensate emerges in the system. In figure [1] we 
present a snapshot of the condensate obtained from numerical simulations for N = 1000 
sites, J = U = 1 and p = 10. The critical density is p c = 0.2397, thus the system is 
deeply in the condensed phase. One sees that the condensate extends over many sites, in 
contrast to the zero-range process where it is always localized. In [TU] the extension W, 
that is the number of sites occupied by the condensate, was estimated to be proportional 
to \/N for sufficiently large systems. The argument was based on a similarity between 
the envelope of the condensate, and a trajectory of a random walker. Here we will 
employ another approach, which will allow us to determine not only the extension, but 
also the average shape of the condensate, its fluctuations and the distribution of particles 
7r(m) in the condensed phase. 

Let us assume that p > p c and consider the weight P n of a configuration, where 
n sites share the surplus of particles M' = N(p — p c ) and where N — n sites form a 
"background" with occupation numbers of order p c . For n = 1, namely for only one site 
carrying the condensate, Pi is approximately 

Pi = NZ N ^(1)K 2 (M' - Pc )p{M' - Pc)- (28) 

The first factor iV stands for iV possibilities of choosing the site occupied by the 
condensate. The second term accounts for iV — 1 "background" sites, it is just the 
partition function ( fl9l) for N — 1 sites. The third term stands for two "domain walls" 
at the borders of the condensate. The fourth term accounts for the on-site potential. 
Since p(m) ^ 1 only for empty sites, the fourth term is simply equal to one. Hence 
the function p(m) influences only Z^-i via the maximal eigenvalue A max . Using the 
definition (fT2|) of K(x) and neglecting p c since it is small in comparison to M' — > oo, 
we can write: 

P l = Nexp[(N - 1) In A max - 2JM'\. (29) 

For n = 2, the condensate may occupy either two neighboring sites, or two condensates 
exist, separated by some number of background sites. The statistical weights for these 




200 400 600 800 1000 300 400 500 600 700 

i i 



Figure 1. Left: a snapshot of the condensate for N = 1000 and p = 10 obtained in a 
MC simulation. Right: the shape of the condensate averaged over many realizations, 
each centered at i = 500. 



two situations are, respectively, 

jingle ~ NZ N _ 2 (l)K(xM')K(xM' - (1 - x)M')K{{\ - x)M'), (30) 

^.double ~ N 2 Z N _ 2 {l)K 2 {xM')K 2 {{\ - x)M'), (31) 

where x is the fraction of particles occupying the first site, and as before we dropped p c 
in the arguments. These formulas can be rewritten as 

^single =P 1 e-*WhJ*'(l-|2x-l|) j (32) 

^.double = Pie" lnW+ln7V , (33) 

so that for any < x < 1, the two-site single condensate with probability P2, single 
dominates in the thermodynamic limit over i) the one-site condensate Pi and ii) the 
case P2, double when the two sites carrying the condensate are separated. 

Similarly, one can show that for n — 3, 4, . . ., P n grows with n. The first observation 
is that the probability is higher when all n sites keep together. For such a case, assuming 
that the ith site takes a portion Xi of the condensate, we obtain 



P n oc exp (-n In A max - JM' \x i+1 - Xi\ + xi + x^j \ , (34) 

where the factor independent of n is dropped, and x± + . . . + x?? — 1. The weight P n 
has a maximum if 

n-l 

f(x) = ^2\x i+1 - Xi\ + Xi + x n (35) 
i=i 

takes the minimum value. To find this minimum, let us observe that f(x) can be 
always lowered if we order x±, . . . , x n so that x\ and x n are minimal, and both sequences 
x±,x 2 , ■ ■ ■ and x n , x„_i, . . . are increasing. Then, x^ has a maximum at some i = i . We 
have: 

io — 1 n— 1 

f(x) = X! + ^2(x i+1 - x^ + - x i+1 ) +x n = 2x io , (36) 

i=l i=io 
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and, since all Xj's must sum up to one, the smallest value is Xj = 1/n. Therefore, 
P n ~ exp(— n In A max — 2JM'/n) clearly grows with n as long as n C \/N and hence 
we see that the condensate must be spread over many sites. One could argue that 
we did not take into account the entropic contribution, namely on how many ways we 
can choose the set of X\, . . . ,x n . But the entropic contribution cannot be larger than 
~ iV n = e nlnN , so as long as n stays finite, it only changes the constant multiplying 
n in (|34p . But if n grows with N as ~ \ N, the two terms ~ n and ~ M'/n become 
comparable. This could suggest that if they both were negative, P n would take its 
maximum at n ~ and the spatial extension of the condensate would be of order 
~ yN. However, the "naive" entropic term n In N dominates over the — n In A max term 
for very large N and therefore configurations with growing n seem to be more probable. 
In other words, this argument would imply that the spatial extension should be ~ N, 
and there would be no condensation at all. The problem is caused by the fact that we 
overestimated the entropic contribution. In the next section we will show that a more 
accurate calculation yields that the entropy ~ n and the extension is indeed ~ y/N. 

3. 3. Shape of the condensate 

The interesting feature of the condensation discussed so far is that the condensate 
is extended over many sites. Apart from the question concerning the scaling of the 
extension with the system size, to which we will return soon, another interesting question 
concerns the shape of the condensate. By the shape we understand here the condensate's 
envelope, that is we ask about the average occupation numbers (m^) as a function of 
the index i. A numerical recipe for measuring the envelope in a Monte Carlo (MC) 
simulation is straightforward: one has to accumulate the shape over many realizations, 
each time shifting the condensation peak to a common origin. The only difficulty is to 
find the center of the condensate. This can be done in many ways which are equivalent 
in the thermodynamic limit. For any finite size, however, we found the following method 
to be the best. First, one searches for the maximal occupation number, let it be the one 
at site i max . Then, starting from i max , one searches for the left ii e f t and the right z r i g ht 
border of the condensate where the occupation numbers drop to the background level 
p c . The center of the condensate is then defined as i = (i\ c { t + Vi g ht)/2. This method 
preserves sharp borders of the condensate seen in figure [U left, quite well, as can be 
seen in the same figure, right, where we show the envelope obtained as an average over 
20000 MC samples. Since steady-state properties of the model do not depend on the 
fact whether it is in equilibrium or not, to simulate the system we treated it as if it 
were in equilibrium with the probability of a microstate (fT3l) . The model can then be 
simulated using the standard MC method by picking up a random particle and moving it 
to a randomly chosen site, accepting or rejecting the move with Metropolis probability. 
This leads to different dynamic properties than the original model has, but considerably 
speeds up the convergence towards the steady state in the condensed phase. 

A careful inspection of figure [1] shows that the envelope is neither circular (which 
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would be true if particles behaved as in a droplet of water due to some surface tension) 
nor Gaussian (the most frequent distribution so that it may be naively expected to 
be observed), but has some parabolic-like shape. A simple argument shows that the 
shape cannot be determined from "energetic" considerations alone, which do not take 
the entropy into account. Imagine that the condensate has a single maximum and 
falls monotonically on both sides. Then the "energy" which is the logarithm of the 
steady-state probability ([TBI and reads approximately — W In A max — J J2i \ m i+i ~~ m i I — 
— W In A max — 2 JH, depends only on the height H and the width W, but not on the 
shape. This means that the principle of energy minimization alone does not allow us to 
estimate the shape. 

Therefore, to find the envelope of the condensate, one has to take the average over 
all possible shapes, weighted by P(m). Below we will show how to do this. Since the 
whole reasoning is valid for more general weights ( TTTT) . as long as K(x) falls off faster 
than any power law, and p(m) = 1 for m > m max , we will keep the discussion quite 
general. At the end, we will present formulas for the special case (TI2"I) . 

Let us remark first that since the grand-canonical partition function ffl5|) does 
not exist in the condensed state, a naive approach via evaluating ( fl5|) for the fugacity 
z determined by (fl6l) does not work. One should in principle work with the canonical 
ensemble, but direct calculation of ffl4|) is very hard. This is because in (fH|) we have p(m) 
which depends on occupation numbers, as well as K(\m — n\) which depends on their 
differences. This prevents us from decoupling terms for different sites and to perform 
the sum over {m^} directly. One can observe, however, that in the condensed state the 
system may be split into two parts: a condensate with W sites, and a background with 
N — W sites. The background can be treated as being at the critical point, because 
the condensate absorbs all the surplus of balls. Summing over all configurations in 
the background we shall get the grand-canonical partition function for the system with 
N — W sites and p c balls per site on average, which reads Z^-wiX) = A^ x w , as follows 
from (Tl9|) . From the formulas 

z dZ N (z) 



(mi) 



NZ N {z) dz 



=i 



(mi(mi - 1)) 



z 2 d 2 Z N (z) 



(37) 



(38) 



NZ N (z) dz 2 

applied to the background we obtain that the variance (mf) — (mi) 2 is constant. Thus 
the total number of particles in the background fluctuates as ~ \/N — W . This in 
turn means that M', the number of particles in the condensate, also fluctuates as 
~ \/ N — W M' . We shall therefore neglect fluctuations of the condensate's mass, and 
assume that it takes always M' particles. The width and the height of the condensate can 
still fluctuate. Moreover, the condensate, from its definition, has occupation numbers 
growing with N. Since we assumed p(m) = 1 for large m, we can drop p(m) in the 
condensate. On the borders, however, « p c and since this is very small in comparison 
with other m/s, we can set it to zero. The statistical weight of the condensate extended 
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to W sites is then 



w+i 



' w 



M' 



(39) 



Z C (W) =J2H K (\ mk ~ mk ~^ 6 

{m k } k=l \_k=l 

where we assumed mo = mw+i = at the borders. Here the function 5[n] denotes the 
Kronecker delta 5 n fl and will be used from now on for notational convenience instead 
of 5 n fl. The total weight of having a system composed of the condensate and the 
background will be assumed as a product of probabilities: 

P(W) = Z N _ W (1)Z C (W) = exp {{N - W) In A max + \nZ c {W)) . (40) 

To find the extension W, one has to find the maximum of (1401) with respect to W. 

We will focus now on calculating Z C (W) defined in ( 1391) . Due to the Kronecker 
delta constraint and the factors K(\rrik — mk-i\) which suppress large differences in 
neighboring occupation numbers, we expect that the majority of occupation numbers 
will be much greater than zero. Therefore, we assume that the summation over {mi} can 
be extended to negative values without changing Z C (W) too much. We will show later 
by calculating average occupation numbers and their variances that this assumption is 
fully justified, that is in fact almost all > in most probable configurations. We can 
now rewrite Z C (W) as 



Z C {W) 



E 

=— t 

w- 

n 



dw+i = — °° 



W+1 



kd k - M' 



k=l 



5 



'W+1 

zZ dk 

.k=l 



(41) 



d\=— oo 
W+1 

X 

k=l 

where = — m^-i- The first Kronecker delta in the above formula gives the 
conservation of particles. The second delta reflects fixed boundary conditions according 
to m = m w+ i = which gives di + . . . + dw+i = 0. Let us introduce an auxiliary 
function 



G(W, u) 



8 



d\=—oo 
W 

X 
k=l 



IV 



^2 kd k 



M' 



k=l 



5 



w 



.k=l 



\\K{\d k \)e d ^\ 



(42) 



where u = {ui, . . . , %} are auxiliary variables. We will use them later, let us now only 
observe that Z C (W) = G(W + 1,0) and hence G(W, u) is a generating function for the 
moments of dk- 



G(W,u) 



-G(W,u) 



u=0 



Replacing both delta functions in (142|) by their integral representations: 

" i7T+t dz „. 



5[x] 



27T+e 



2ni 



(43) 



(44) 
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and performing the sum over {dk}, we obtain: 



G(W, u) 



t7T+€2 



dv 



where 



F(z, v, u) = -M'z + ^ In K(u k +v-kz). 
The function K(x) is defined as 

oo 

k(x) = j2 k mv 

d=—oo 



fc=l 



(45) 



(46) 



(47) 



The function G(W, u) can be evaluated in the saddle point: G(W, u) ~ e p { z ^^ u ) ^ 
with z = z(u),v = v(u) being solutions to the saddle-point equation d z F(z,v,u) = 
d v F(z,v,u) = 0. Let us assume for a moment that we have solved this equation and 
have determined both z(u) and v(u). With help of (j430 . we can then write: 

" d 



(4) = 

var(4) 



F(z(u),v(u),u) 



du k 

- (4) - (4) 



u=0 



d«2 



F(z(u),v(u),u) 



J u=0 



cov(dj,d k ) = (djd k ) - (dj) (4) 



d d 

duj duk 



F(z(u),v(u),u) 



(48) 
(49) 
(50) 



M=0 



and similarly for higher-order correlation functions. We stress that in general all 
derivatives have to be taken for z(u),v(u) being functions of u, and only at the end 
one can set u = 0. But for the first moment we have 



(4) 



OF dz OF dv OF 
+ "n 7^ + 



(51) 



u=0 



dz duk dv duk duk 
where the partial derivatives of F are taken in the saddle point and thus are zero, except 
of the last one. We thus obtain: 

K'(v - kz) 



(d k ) 



(52) 



K{v — kz) 

The first observation that eliminates z without the need of solving the saddle-point 
equation is that the averaged peak must be symmetric around its center. This is true 
even if the system is not in equilibrium, that is when a current of particles flows through 
the system for r ^ 1/2. The symmetry of the shape implies that (m^) = (mw-k), thus 
(4) = — (dw-k) and hence z = (2/W)v. We obtain: 

K'(v(l - §)) 



(4) 



(53) 



In addition, because the total number of particles in the condensate equals M', we have 
for large systems: 



X><4> 



fc=0 



w k' (v (1 - f )) 
„ k K Hl-f))* 



-M'. 



(54) 
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The last equation, after changing variables, reduces to 



w = v 



xK'(x) 



dx 



-1/2 



(55) 



where neither M' nor W appear alone and where we defined the reduced extension 
w = W/ V M'. Equation (1515]) fixes w as a function of v. We can now write 

2v 



In Z C (W) F 



v,0 



(56) 



and combining equations ([3611 . (H6l) and (jlQl) we obtain that the logarithm of the weight 
of a condensate extended over W sites is given by 

2v 



\nP(W(v)) = N\n\ m 

w(v) 



M' 



-w(v) In A r 



W[V\ 



lnK(x)dx 



(57) 



Since does not appear in the above formula explicitly, but only through the relation 
W = w(v)\/M', the extension W is now determined by maximizing fl57|) with respect 
to v. Hence, W must be proportional to yM' ~ \/N. Thus, we have proved that the 
extension grows with the square root of the system size. 

One can go one step further and calculate v and w, and thus the exact, numeric 
value of the extension W. The probability (IBTj) assumes a maximal value, when its 
derivative with respect to v is zero. Inserting fl55l) into fl57j) and taking the derivative 
d/dv we obtain: 



In A max - In K(v, 



1 nif'W 



0. 



(58) 



v 2 K(t») xK ' (x U x 

One can show (see Appendix A) that the second term in the square brackets is never 
zero. Only the first term counts and gives: 

v = ^ _1 (A max ). (59) 

Having v and w, it is now easy to obtain the shape of the condensate: 



\fc=l / fc=l 
where (<ifc) is given by fl52|) . If we define new variables: 

2n 



= (m n ) /VM>, 



t 



- 1, 



so that t G [—1,1], we can write the formula for the envelope h(t) as follows: 

, . . w , K(v) 
h(t) = — ln^ v 



2v K(vt) 



(60) 



(61) 



(62) 




Figure 2. Left: Comparison between the envelope of the condensate from (|64|) and 
rescaled MC data for the weights (U]) with J = U = 1 for N = 1000, p = M/N = 1 
and 3 (circles, squares) and N — 4000, p — 1,3 (diamonds, triangles). Right: the 
envelope (|64|) (solid line) compared to the parabola (|66|) (dotted line). 



Thus the shape in reduced variables t, h(t) becomes independent of the system size and 
the density of particles. The parameters w, v depend on the weights K(x) and p(m) 
only. 

Let us now go back to the choice (Tl2i) of the weight functions. For this case, one 
obtains 

K(x) = sinhJ (g 3 N 
cosh J — cosh x ' 

and v, calculated from (1591) . assumes a very simple form: v = J — Jq. Inserting v 
obtained from the above equation into (I5"5T) . we can calculate w numerically. The shape 
is given by (1621) and reads 

,.. w, ( cosh J — cosh vt \ . . s 

h(t) = —ln( — — , (64) 



2v \ cosh J — cosh v 

with w, v being functions of the parameters J, U . For instance, for J = U = 1 one 
obtains v = 0.5413 and w = 2.2005. In figure [2] we show a comparison between (1641) 
and MC simulations for different system sizes. One sees that the curves obtained for 
different N, p and plotted in the rescaled variables t, h agree with the theoretical line. 
Small deviations in tails can be attributed to finite-size effects and will be discussed in 
section 3.4, but one sees already in figure [2] that slopes of the envelope become more 
steep for increasing system size N. 

One may ask now the following question. The shape of h(t) is determined by (14T]) . 
which is just the partition function for a Brownian excursion: a ID random walk with 
steps drawn from K(x) starts at m = and ends also at m = after W steps. The only 
difference in comparison to ordinary excursions is that the area under the trajectory 
is fixed to some value M' . If this area is not fixed, it fluctuates around ~ W 3 ^ 2 , the 
result known for a "free" excursion [T7j. Assuming M' cx W 3 ^ 2 in our formulas we have 
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w = W/ \/M' oc W 1 / 4 and from (j55p we obtain that 

V = 1^L^ . ( 65) 
K"(0)w 2 

Thus, by applying Taylor expansion to (1521) one obtains a universal curve, independent 
ofK(x): 

h{t) = (6/4w)(l -t 2 ). (66) 

It must be stressed, however, that in our case M' scales always as W 2 , so the above 
argument does not hold, and the envelope is not universal. In figure [2j right, we see that 
despite a lack of universality in general, the exact form (IMj) does resemble the parabola 
from (1661) very much, at least in a certain range of J,U. 

3.4- Fluctuations 

In this section we shall consider only the special choice ( 1121) . Assuming that the width 
W of the condensate is fixed, the fluctuations around the average envelope are measured 
by: 

n n 

var(m n V = ^^((dkdj) - (d k ) (dj)) 

fc=l j=l 



n 



var(4) + 5^52 cov(4, dj), (67) 



fc=i fe=i jyfc 



where the subscript tells us that W is assumed to be constant. We need to calculate 



8 2 F 8 2 F ( dz\ 2 d 2 F ( dv ^ 2 

T 

k 



dul dz 2 \duk) dv 2 \duk 



and 



^ d 2 F dz ^ 9 2 F dv ^ d 2 F dz dv 
dzduk duk dvduk du k dzdv du k du^ 



d 2 F dz dz d 2 F dv dv d 2 F ( dv dz dv dz 

CQV(4, dj) = -TrTTT- 7T- + -77^77- + 7^7 [^1^7 + 



<9z 2 <9t> 2 c^Mfc <9wj <9z<% \ duk duj duj duk 

d 2 F dz d 2 F dz 
+ n n n + 



dzduk duj dzduj duk 
+ d 2 F dv + d 2 F dv ^ 
dvduk duj dvduj duk ' 

where all partial derivatives of F(z, v, u) defined in (H61) are taken in the saddle point 
with u — 0. Other derivatives which in principle may appear in the above formula are 
zero and have been dropped. Since this is rather technical, we postpone the derivation 
to Appendix B. Here we shall only present the final result: 
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where the size-independent variance cr 2 (y) is given by 

i 2 (y) il(y) 



a\y) = I Q (y) - 



with 



I m {t) 



y" 



cosh J cosh {vy) — 1 



dy 



(71) 



(72) 



. i (cosh J — cosh (vy)) 2 

derived in Appendix B. The formula (ITUj) would give the fluctuations of the envelope 
if the width W was kept fixed. In reality, the width also fluctuates around the mean 
W = M' . To take this into account we have to sum over all possible widths, each 
having the probability P(W) given by (157|) . 



var(m n ) = J>W ( m Dw ~ ( < m ™>' 
w V w 

= ^P(W) [vai(m n ) w + ((m n ) w - (m n )) 2 ] 



(73) 
(74) 



Here (m n ) w means the envelope for given W, while (m n ) means the envelope averaged 
over W. Since the probability P(W) is concentrated around W = wy/ M' for large 
systems, the latter can be approximated by V M'h(t), where t = 2n/(w\/ M') — 1. 
Also the variance vai(m n )w summed with weights P(W) is well approximated by 
7;W\/M'a 2 (t). Then, from the above equation it follows that 



var(m(t)) 



'M' 



^a 2 (t) + y/M' I P(W) [hit 



wVM' \ wVM' 



2 w ./ \ \ ir / ir 

-h(t)) 2 dW] , (75) 

so there is an additional contribution to the variance from fluctuations of the 
condensate's width. The probability P(W) can be well approximated by a Gaussian 
distribution with the mean wVM* and variance var(W) = V M's 2 with 



■K'(v) - K(v) 



s =w— = . 76 

vK'(v) 

This formula comes from expanding the probability P(W(v)) around the minimum for 
v = Vo, and is true for any K(x). For our special case ( fl~2]) it reads 

w 2 /4 + (cosh(t>) — cosh (J))/ sinh(f) 



w- 



(77) 



Then, equation fj75j) can be rewritten as 
var(m(t)) w vM 



w 2/ , «)M' 3/4 
77 * 2 (t) + 



X 



V2tts 



(7f 




Figure 3. Left: The variance var(m n ) obtained from MC simulations for N = 
4000, p — 3 (circles), compared to (|78j) evaluated numerically (thin line). Right: the 
averaged envelope (h(t)) from taking into account small fluctuations of M'. 

In figure El left, we compare the result (175]) and the one obtained by means 
of MC simulations. To a good approximation, the variance var(m n ) grows linearly 
with W ~ \J M'. This means that in the thermodynamic limit the number of 
occupation numbers being less than zero becomes negligible, because typically m n = 
(m n ) ± y/var(m n ) ~ W ± yW, so that extending the summation in (j41~|) to negative 
values was fully allowed because it leads to self-consistent results. At the center of the 
condensate, for t = 0, the formula (ITS]) simplifies to 

var(m(0)) » VM 7 -a 2 (0) + -^/i 2 (0) . (79) 
\ 2 k;^ / 

In the next section we shall see that this formula is very important for the problem of 
the distribution of particles. 

As a by-product we obtained the variance V M's 2 of the distribution of the width 
W. This allows us to calculate the envelope h(t) averaged over fluctuations of W: 

wJgJ'^H^ (80) 

In figure M, right, we show the curve obtained from the above equation. It fits to the 
MC data points even better than equation (IMl) . especially in the tails. This indicates 
that smooth tails are mainly due to the fluctuations of the condensate's width. This 
effect clearly vanishes with increasing system size, i.e., is a typical finite-size effect. 

It turns out that not only (m n ) and var(m ra ) but also the distribution of m n for 
any fixed n can be predicted. Since m n is a sum of many (almost independent) random 
variables d^s, each of them having finite variance, the distribution tends to a Gaussian, 
according to the Central Limit Theorem, provided that l<n< W, that is we are not 
too close to the borders. In the thermodynamic limit we can safely assume that all m n 
have Gaussian distributions. 
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3.5. Distribution of particles 



Let us denote by 7Tb (m), 7r c (m) the probability of finding m particles at a background 
site and a condensate site, respectively. The distribution 7Tb (m) oc 0^, as follows from 
the partition function (TT5I) at the critical point. In particular, for the previous weights 

2m 



one obtains 7Tb (m) oc 



oUS m ,0 



The distribution in the condensate, ix r 



771) , IS 



given by summing Gaussian distributions with mean (m n ) and variance var(m n ) over 



all n 



W condensate sites: 



7r r m oc 



w 

E 



' y / 27rvar(m n ) 



exp 



(m„)) 2 



2var(m n ) 



which gives 



7r c (m) oc 
By changing variables t 



dt 



= exp 

o A/27rvar(m(t)) 

y = h(t) one obtains 



(m- VWh^f 
2var(m(t)) 



52) 



7r c (m) oc 



h(0) 



(O'(y) 



exp 



(to-v 7 !? 7 ^) 2 
2var(m(/i -1 (3/))) 



■dj/. 



83) 



A/27rvar(m(/i 1 (y))) 

Hence the probability 7r c (m) may be viewed as a convolution of (h^ 1 )'(y) with a 
Gaussian distribution which smears the thermodynamic-limit distribution 7r c (m) ~ 
{h~ l )'(m/ \J M'\ Since the variance grows like y/M' , the smearing acts on distances 
~ M' 1 / 4 and influences the profile of 7r c (m) only at y ~ /i(0), because (h~ 1 )'(y) is 
narrow only at y = /i(0). We can therefore assume var(m(t)) = var(m(0)) because it 
affects the distribution only at t = 0, and write 



7r c (m) oc 



dt 



: exp 



(m - \fMh(t)Y 
2var(m(0)) 



4) 



'o A/27rvar(m(0)) 

with var(m(0)) given by (!79|) . Since the main contribution to 7r c (m) comes from the flat 
region in h(t), the condensate's peak tends to a Gaussian for large systems. 

The probability 7r(m) of finding m particles at any site is additively composed 
according to 



7i{m) 



N 



-W W 



55) 



In figure H] we show a comparison between 7r(m) obtained from (!8~5j) after inserting 7Tb (m) 
and 7r c (m) normalized to probabilities, with MC data. The agreement is very good. 



4. The case when p(m) ^ 1 for all m 

As long as K(x) falls off faster than any power-law, and p(m) = const for m larger than 
some m max , the method from the previous section directly applies. This means that 
one can calculate the envelope of the condensate, its fluctuations and the distribution 
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Figure 4. Left: the distribution of particles 7r(m) obtained from (l85|) (thin line) and 
from MC simulations for N — 8000, p ~ 3 (circles). Right: the same on a log-log scale. 



of particles. Since p(m) can be multiplied by any factor without changing "physical" 
quantities, the condition on p{m) can be assumed to be p{ra > m max ) = 1. 

However, if p(m) does not tend to one for large m, the partition function cannot 
be approximated by (1391) . An important example is 



K{x) ~ e~ aW , p(m) ~ e- bm \ (86) 

for a,b, /3,7 > 0, with both K(x) and p(m) decaying fast with x and m, respectively. 
In order to have the condensation, the grand-canonical partition function must have a 
finite radius of convergence. This means that the function 

9(™, n)z m = e~ bnl e~ a \ m - n ^- bm \ m , (87) 



must diverge for z larger than some critical z c . It is easy to check that z c is finite for 
7 < 1 and there is condensation above some critical density. For 7 > 1, the radius of 
convergence is infinite, so there is no condensation for 7 > 1. The value of (3 is not 
important here and can be arbitrarily large. 

Next we have to determine for which values of (3, 7 the condensate is extended. Let 
us calculate the contribution from the condensate having ~ iV particles at a single site. 
This is the same approach as in section 13.21 

P, W Nc^p^K^N) W Nc N-l e -bN^-2aN^ (88) 

As before, N stands for iV possible condensate's positions, c N_1 with some constant c 
accounts for background sites, p(N) is the on-site potential at the condensed site and 
K 2 (N) stands for the two domain walls. The contribution from the condensate residing 
on two adjacent sites, having respectively N/2 + e and N/2 — e particles, reads 

P 2 w Nc N ~ 2 K(N/2 + e)K(N/2 - e)K(2e)p(N/2 + e)p(N/2 - e) (89) 



Nc 



N-2 e -a(N/2+e)' 3 -a(N/2-e)' 3 -a\2e\l 3 ~b(N/2+e)^-b(N/2-€yi 



(90) 
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The contribution is large only for e ~ 0. We can then write: 

In P 1 /P 2 « - In c - 2o(l - 2~ P )N P - b(l - 2 1 " 7 )iV 7 . (91) 
Since 7 < 1, 1 — 2 1-7 < 0. On the other hand, 1 — 2 _/3 > and hence 

In Pi /.P 2 - -N p + (92) 

where the proportionality factors have been skipped. For 7 > (3 we thus obtain 
P1/P2 — ► 00, so the condensate must be located on a single site. For 7 < /3, however, 
the contribution from two sites is larger than from a single one. Similarly, one can 
show that the contribution grows for 3,4,5,... sites. Therefore the condensate must be 
extended, but we still do not know, for how much. 



4-1- Fixed-envelope approximation 

We have learned that fluctuations of the occupation numbers can be neglected if 
p(m) — > 1. This is also the case here, because the fluctuations are determined by 
the variance of K(x), which is again finite. We shall thus assume that, for sufficiently 
large systems, the envelope of the condensate is essentially fixed, modulo some small 
fluctuations around it. Then, the logarithm of the weight of the condensate extended 
to W sites is 

In P(W) « - Win A max + Ws + J2^K ((m k+1 - m k )) 

k 

+ £}lnp«m fe )). ( 93 ) 

k 

Here s is some "entropic" factor due to small fluctuations of m^'s. Values of (m*.) are 
equal to Hh{2k/W — 1), where H = M'/W is the height of the condensate and h(t) is 
some fixed "envelope" having the same meaning as in section [373]. but its functional form 
may be unknown. The differential term (m^ +1 — m k ) is more difficult. It can behave in 
two distinct ways. First, if we assume that h(t) is a "smooth" (differentiable) function, 
the difference can be rewritten as ^rh'(t) with t = 2k/W — 1. A second possibility 
is that h(t) has a rectangular shape, h'(t) does not exists at t — ±1. Thus, we may 
have in principle two extended condensates: the one whose shape resembles that in the 
previous sections, with smooth h(t), which we will call "smooth", and a "rectangular" 
one where h(t) looks like a step function in the limit of large N. We can write 



lnP mM (W)^W 



+ I lnp{Hh{t))dt 
Jo 



(94) 



lnP Tect (W) w Wc + 2\nK (h(l)^\ + W J^ lnp(Hh(t))dt, (95) 

where c = s — lnA max is some unknown constant. To find the extension W, we must 
find the larger of the two maximum values of the P(W / )'s. Assuming now the weights 
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with < 7 < 1, (3 > 7, we have: 
In P smooth (W)^W 



o 



7 /•! 



-M^J / IM*)l 7 d* 



(96) 

P / Jf \ 7 rl 



fM'\ p f m'\^ r 

lnP rect (W)~Wc-a( w ) 2\h(l)f-bW[ w ) J MtWdt. (97) 

The integrals are just some constants. The value c must be smaller than zero, because 
otherwise P(W) has no maximum which is in contradiction to the fact that the 
condensate is extended. To estimate the extension, we can drop constants (remembering 
that they are all positive, except c) and write 

fM'V fM'V 
lnP smoot hW~ -W-Wl—j -W( — \ , (98) 

fM'V fM'V 
\nP TCCt (W) --W-f^) - W (W) ■ (99) 

Let us consider P smoo th(W) first. If (5 > 1/2, it has a maximum for d\nP smoot ^/dW = 
which reads 

fM'V fM'V 
(1-2/?) (— -) + — =-1, (100) 



W 2 J \w 
and hence 

W ~ M'V > -rVW-y) } ( 101 ) 

because the exponents of M' in both terms must be equal. We also obtain the maximal 
value: 

In PsmoothW ~ -MW-tV( 2 ^ . (102) 

If (3 < 1/2, the above equation has no solution, and In P smooth (H / ) takes its maximal 
value at W ~ 1, which would mean "delta-like". This is impossible, so this solution is 
obviously not valid for j3 < 1/2, where the condensate must therefore be rectangular 
and described by (|9~9~I) . which has the maximum for 

pM'^W-P- 1 - (1 - 7)M'W- 7 = 1. (103) 

From the above equation we obtain the extension of the rectangular condensate 

W ~ M'V-^W-i+V , (104) 

as well as 

In P Iect (W) ~ -M ,/3/(/3 " 7+1) . (105) 

If /3 > 1/2, both equations (|98l) and (1991) are valid. But for M' —+ oo, it turns out that 
lnP rect (iy) > In P smooth (W / ) for 1/2 < (3 < 1. Therefore, collecting everything together 
we see that the rectangular-shaped condensate dominates for f3 < 1, and the smooth 
condensate dominates for (3 > 1. In figure [5] we show a phase diagram which summarizes 
all the above results. One sees there that the extension W ~ M' a can be tuned to any 
< a < 1/2, but cannot be larger than 1/2. 
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Figure 5. Phase diagram for K(x) ~ e - ' 31 and p(m) ~ e - ™ 7 . Values of the exponent 
a in the extension W ~ M /a , a rcc t = {fi — l)/{fi ~ 7 + 1) f° r the rectangular and 
"smooth = {fi — 7)/(2j9 — 7) for the smooth condensate, are represented by the colour 
(gray) code. Dotted lines mark a — 0.05, 0.1, . . . , 0.45. 



5. The case when K(x) has heavy tails 

In this section we consider K(x) decaying as a power-law: K(x) ~ x~ v . For this case, 
the analytic method from Section 3 cannot be applied even if pirn) is constant above 
some m max . The reason is that the series (|4"Tj) does not converge for K[x) decaying more 
slowly than exponentially. On the other hand, the assumption that fluctuations vanish 
in the thermodynamical limit, a crucial step in the fixed-envelope approximation, is also 
not well justified a priori, because K(x) may have infinite variance. Fortunately, it will 
turn out to be quite easy to show that the condensate is never extended over more than 
some finite number of sites, so there will be no need to use any of the two methods 
mentioned. 



5.1. Double power law: K(x) ~ x v and pirn) ~ m b 

Let us first discuss the criterion for condensation. For b > 0, the eigenvector m has the 
following form for large m: 

<P m ~ m- fe / 2 ^, (106) 

which can be checked by inserting it into (!23j) . The critical density 

P c~Y. m ~ b ~ 2v+l ( 107 ) 

rn 

is finite when v > 1 — 6/2. The above formula is valid also for v < 0, that is 
when configurations with large differences between neighboring occupation numbers 
are favored. For v = 0, that is when K(x) ~ const, we recover the ZRP case: the 
critical density is finite for b > 2. If b < 0, the solution (11061) is not valid, but then m 
has to grow with m. The critical density is therefore infinite for b < 0. 

To estimate the extension, let us first calculate the contribution from the condensate 
having ~ iV particles at a single site. Following the same method as before one obtains 

Pi w Nc^p^K^N) w c N - l N l - b - 2v . (108) 
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The contribution from the condensate residing on two adjacent sites, taking respectively 
x and 1 — x fractions of particles, is: 

P 2 » Nc N - 2 p(xN)f ((1 - x)N) K(xN)K ((1 - x)N) K (|1 — 2ac|iV) .(109) 

For fixed x, in the limit iV — > oo, we obtain 

P 2 w c^iV 1 - 26 - 3 "^^ - x))- 6 -"|l - 2x1-". (110) 

Clearly P 2 ~ PiN~ b ~ u and since —b — v < —6/2 — 1, it vanishes in comparison to Pi 
in the thermodynamic limit. Even if we assume that the difference e in occupation 
numbers of both sites is small (x ~ 1/2), we have 

P 2 w c^iV 1 - 26 - 21 '^), (111) 

which is still by a factor ~ N~~ b smaller than Pi. Performing also the sum over e gives 
an additional factor not larger than N 1_v . Thus P 2 ~ P L N 1 - b - v < PiN~ b / 2 and the 
ratio P 2 /Pi tends always to zero. The conclusion is therefore that the condensate, if it 
exists, is always localized at a single site for any b, u, precisely as for the ZRP. 

This localization is, however, not easy to observe in numerical simulations. If the 
system is initially prepared in a state with a uniform distribution of particles and the 
simulation starts, after even a very long time one usually sees the condensate occupying 
more than one site. This is caused by a small probability to decay from the extended 
to a single-site condensate. Once the condensate becomes localized on a few sites, there 
are large energetic barriers between condensates of decreasing sizes. On the other hand, 
if one starts the simulation from the state with only one site occupied by all particles, 
one observes that this state does not decay to an extended condensate in the course of 
the simulation. 



5.2. Mixed decay law: K(x) ~ x " and p(m) 



MS, 



m,0 



Let us finally briefly discuss the case where K(x) ~ x~ u follows a power law as in the 
previous subsection but p(m) = e USm >° is of short-range nature as in (fl2"|) . Again, we 
shall first estimate the probabilities Pi, P 2 of having the condensate on one or two sites: 

Pi ~ Nc N ~ l N~ 2u , (112) 
P 2 ~ Nc N ~ 2 N- 2v A v K(e). (113) 

One sees that both terms depends on A" in the same way and none becomes negligible 
for A" — ► oo. Also higher P n share the same feature: 

-2 



Nc 



N-n 



(114) 



This suggests that the extension is determined by the ratio K{x)/c (or more 
precisely: by the inverse of its logarithm) and is constant if this ratio is smaller than 
one. If K(x) = oo, that is for v < 1, or when ^ x K(x) /c> 1, P n grows with n — > oo 
and the condensation is impossible. This means that the prefactor in K(x) ~ x~ v does 



Mass condensation in one dimension with pair-factorized steady states 



26 



matter — consider for instance K{x) = 1/(1 + x/a) u which behaves clS Qj^ ' X v for large 
x. When a is small, Yl x K{ x ) is large and there is no condensation. 

We can get the same prediction by applying naively the fixed-envelope approach, 
developed in section HI and forgetting about problems caused by (perhaps) large 
fluctuations of neighboring occupation numbers. We use equations and ( 1951) with 
K(x) ~ x~ v and pirn) = 0, because p(m) vanishes now in the condensate. We find that 
lnP smooth ~ — M 1 / 2 with W ~ M 1 / 2 and lnP rect ~ —const with W ~ const, so that 
In P sm ooth < In -Prect and the condensate is rectangular and of fixed extension. 

6. Conclusions and Outlook 

We have shown how to tune the shape of the condensate during spontaneous symmetry 
breaking in mass transport models on a ring topology. The shape of the condensate and 
the scaling of its extension with the system size are non-universal. They depend on the 
competition between the ultralocal and local contributions (in occupation number space) 
to the weight factors whose product over all pairs of sites determines the stationary 
state. Analytical predictions were possible even for the condensed phase above the 
critical mass density for which the grand-canonical partition function is no longer 
convergent; a partition function that factorizes over the condensate and the critical 
background turns out to be the appropriate approximation scheme, as the excellent 
agreement with numerical simulations of the shape and the single-site mass distribution 
have demonstrated. 

Some of the results presented in this paper we have already used in [16] for 
predicting the onset of condensation in case of anisotropic hopping in two dimensions, 
where the two-dimensional system could be dimensionally reduced to an effectively one- 
dimensional ZRP. Moreover, in [16] we generalized the topology on which certain classes 
of hopping rates lead to PFSS from a one-dimensional ring to arbitrary connected and 
undirected graphs. From the theoretical point of view it remains challenging to derive 
the phase structure from these known PFSS in higher dimensions that cannot be reduced 
to effectively one-dimensional processes. 

Vice versa, from experimental observations of the shape and the scaling of the 
width with the system size, one may trace back the class of hopping interactions that 
are compatible with the observations. When atoms condense on a crystal surface, they 
can migrate and build extended islands. As experiments on the deposition of clusters [H] 
or fabrication of quantum dots [15] show, the islands can be extended in the direction 
perpendicular to the surface. Currently it is open as to whether hopping rates leading 
to PFSS are able to reproduce the shape and the typical size of the islands of atoms 
obtained in such experiments. 
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Appendix A 
We will show that 

entering (1581) is never equal to zero. From the definition of K(x), 

oo 

K(x) = £ K(\d\)e dx , (A.2) 

d=— oo 

we see that K'(x) + K(0) > K(x) for x > 0. This in turn means that K'(x)/K(x) 
grows with x. Therefore, xK'{x)/K(x) is convex for x > 0. The area under a convex 
function f(x) over a range (0, v) is smaller than the area of a trapezoid whose parallel 
edges are placed at x — 0, v and range from zero to /(a), fib), respectively. This leads 
to the following inequality: 

€M d , < t m, (a.3) 

K{x) 2 K{v) 



and finally gives 



„ 1 vK'iv) 1 /A , 

Q < — -tt = 0, A.4 

v 2 K(v) 

V 7 2 K(v) 



so that Q < for all v > 0. 



Appendix B 

We shall evaluate the formulae fl68|) and (1691) in order to find the variance vax{m n )w 
from ( 1671) . They can be simplified if one observes that by differentiating the saddle point 
equation with respect to one obtains: 

d 2 F dz d 2 F dv d 2 F 
dz 2 dui dzdv dui dzdui ' 

d 2 F dv + d 2 F dz _ d 2 F ^ 
dv 2 dui dzdv dui dvdui 
One can insert equations (IB. II) . (IB.2[) into (16S1) and (1691 . and rewrite (I6T|) as 

^ dui ^ ^ \dzdu k duj dvdu k duj ) ' 

k=l K k=l j=l J J ' 



Mass condensation in one dimension with pair-factorized steady states 



28 



This formula is far less complicated, but still needs derivatives of z, v taken at the saddle 
point with u — 0. Instead of solving the saddle-point equation with directly, and 

differentiating the solution over Ui, we can use the following trick: we will calculate the 
derivatives of F and obtain dz/dui, dv/dui from equations (IB. II) . (IB. 21) . Let us calculate 
first: 

d 2 F 

SSL =c(k) ' (B ' 4) 

d 2 F k 

- -c(fc), (B.5) 



dzduk zq 
d 2 F 1 



., -c(fc), (B.6) 

f) 2 F 1 1 

' ' - 2 k 2 c(k) = ~ 2 -^-(/o(l) + / 2 (1)), (B.7) 



dz 2 z 2 ^ z 2 
d 2 F 1 1 W 



1 i w 



dv 2 v 2 ^ v 2 2 



FP 1 F 1 1 M/ 2 

^ =-— EM*)-- 7 W, (b.9) 



k 



where we defined 



, x cosh Jcosh \v (l - #)1 - 1 

C(k) = ^ ^ j, ( R1 °) 

[cosh J - cosh [y (l - fp))] 
T / x /"* ^ cosh Jcosh(fy) — 1 , 

f -" |s J/ Mj-^> (R11) 

From (1B.1I) and (1B.2I) one sees that it must be also J^- ~ c(i) and Jj^ ~ c(i). The 
proportionality factors in these formulas must be functions of type A + Bi, because the 
only dependence on i in equations (IB.lj) . (IB. 21) . after dividing by c(z), is linear in i. 
Inserting = (A z + B z i)c(i) and JJr = (A„ + -B„i)c(i) into flB~H and flB~2l) and solving 
for A V ,A Z ,B V ,B Z , we obtain 

* -TT^fl-lW), (B.12) 



/ 2 (i)w 2 V w 

2 ' 1-^ + 1)*). (B.13) 



d Uj h{l)W \ 7 (1) W 
Inserting this into (1B.3|) . we finally obtain the variance for fixed width W: 

, s wx/M 7 9 / In \ , 
var(m„V = — - — o 2 — = - 1 , (B.14) 



2 \wVM' 
where we defined a size-independent variance a 2 as follows: 

This is precisely equation (I7TT) . 
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